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We compute the transmission properties of 2-D electromagnetic TM waves that 
are normally incident on a Fabry-Perot structure with mirrors consisting of pho¬ 
tonic crystals. We use a boundary integral formulation with quadratic boundary 
elements and utilize the Ewald representation for the Green’s functions. We trace 
the frequencies of the Fabry-Perot cavity modes traversing the photonic bandgap 
as the cavity length increases and calculate corresponding Q-values. For the case 
of lossy dielectrics, we compare bandgap frequencies and Q-values to experimental 
results obtained by H. Everitt and his group at Duke University. 


1 Introduction 

We model the scattering of electromagnetic waves in 2-D photonic crystal 
lattices consisting of one dielectric material embedded in another in a config¬ 
uration which is spatially periodic in the X and Y directions and homoge¬ 
neous in the Z direction. For these structures, the lattice geometry, dielectric 
contrast and dielectric loss properties of a photonic crystal all influence the 
manner in which electromagnetic waves propagate through the crystal. It is 
well known that 2-D crystals exhibit photonic bandgaps , that is, intervals in 
their frequency spectrum over which there are no waves propagating through 
the lattice (e.g. Joannopoulos et al. 5 , Yablonovitch 13 ). Localized defect modes 
can appear as resonant peaks at bandgap frequencies when local physical de¬ 
fects are introduced into the 2-D infinite lattice. 

In practice, photonic crystals have finite extent which plays a central 
role in the design of optical and electronic devices including filters, lasers 
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Figure 1. Fabry-Perot Structure with Photonic Crystal Mirrors. 


and microwave antennas. Quantitative characterizations of transmission near 
resonant frequencies are needed for optimal design of practical devices. For 
example, the efficiency of a laser cavity (built as a defect in a photonic crystal) 
increases as the transmission peak at the resonance becomes sharper. 

We consider finite extent by truncating the infinite 2-D periodic lattice in 
the X direction so that the dielectric permittivity is constant outside a certain 
interval in X and we model the scattering of plane electromagnetic waves. We 
employ a boundary integral formulation which reduces the problem from 2- 
D to 1-D and also enforces, automatically, radiation conditions through the 
choice of appropriate Green’s functions. The Ewald 4 representation of the 
Green’s function is used to ensure rapid convergence and a further simple 
reduction is applied to facilitate numerical evaluation. Our formulation is valid 
for both lossless and lossy dielectrics, for waves incident at an arbitrary angle 
and is flexible with respect to the geometry of the dielectric configuration. It 
is described in detail in a recent paper. 12 

In this study, we model an experimental setup of H. Everitt and his group 
in the Physics Department at Duke University. We analyze the transmission 
properties of a plane transverse magnetic (TM) wave that is normally incident 
on the geometry of Figure 1 (J — 6 in the figure). The geometry consists of 
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two square lattices of parallel circular rods, truncated in the X direction, 
that are separated by a Fabry-Perot cavity of length d > 0. In particular, we 
characterize defect modes that are of central interest to the experimental work 
in the case of both lossless and lossy dielectric rods. Although we limit the 
present study to circular rods and normal incidence, our approach is applicable 
to a wide variety of 2-D geometries with periodicity in the Y direction. 


2 Governing Equations 
2.1 Maxwell's Equations 

For a time-periodic electric field E(R)e _iw ‘ and a time-periodic magnetic field 
H(R)e _ ’" t , the governing equations in both the interior and exterior of the 
dielectric rods are Maxwell’s equations: 

V x E(R) = —H(R), V x H(R) = -^e(R)E(R) (1) 

c c 

V • H = 0, V -E = 0 (2) 

where R = (X,y), c is the speed of light and e is the relative dielectric 
permittivity of the medium. In the case of E polarization (transverse magnetic 
or TM waves), E = (0,0, VO and H = {H u H 2 ,0) and Maxwell’s equations 
can be reduced to the Helmholtz equation for ip: 

&+f%)1>+^<X,Y)iP = 0 < 3 > 

where the permittivity e(X,Y) is periodic in the variable Y . In addition to 
this equation, we have matching conditions on the rod surfaces (dielectric 
interfaces) which require continuity of ip and Vip. In considering (3), the 
following conditions apply: 

• The field ip in the exterior is the sum of an incident plane wave ip in (X , Y) 
and a scattered wave il>. c {X,Y), that satisfies (3) plus a radiation condi- 
tion at x = ±oo. 

• On the interface of each rod, the following matching conditions are sat¬ 
isfied: 

Ip = V'in 4* V»*c = IPinU 9 n 1p = d n 1p in -f d n 1p» c = AnV'int (4) 
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2.2 Nondimensionalization 

We assume a Y-periodic geometry of period a and plane wave incident rar 
diation of frequency /. We introduce the nondimensional space variables 
x = 2nX/a y y = 2 tt Y/a and the reduced frequency k = /a/c = u)a/(2nc). In 
the new variables, the Helmholtz equation for the field ip is: 

(d? + tf)il> + k 2 e(x,y)il> = 0 (5) 

where we have retained the symbol ip for the field and the permittivity function 
e(x,y) is 27 r-periodic in the variable y. In this paper, we assume that the 
exterior dielectric permittivity is unity, so that: 

/ 1 outside the rods 

e \ x i y) — ^ £ * > i inside the j -th rod 

The quantity £j is thus the dielectric contrast £j = UntjUxt of the jth rod. 
For this study, we restrict our analysis to normally incident fields which have 
the form: 

1 \>in(x,y) = e ikx ( 7 ) 

We allow for the consideration of both lossless materials ( k real) and lossy 
materials in which k can have a small positive imaginary part. We require 
that the scattered field ip BC also satisfy the pseudo-periodic condition: 

ip»c(x y y 4- 2tt) = ip, c (x , y) (8) 

As a consequence of this periodicity, we can restrict our attention to a single 
interval of length 27r in the y-direction. 

3 Boundary Integral Equations and Green’s Functions 

3.1 The System of Integral Equations on the Boundary 

We denote the jth dielectric interface by ODj and the union of all interfaces by 
9D = Vj„idDj. The governing equation (3) can be reduced to the following 
system of integral equations on the boundary which must be satisfied for each 
f G 9D : 12 


y. c (v)~ [ *£4 + f G{r-r)d n Mr)d,(r) = 0 
2 JOD JoD 


(9) 


+ / 9 *£- r k . c (r)ds{r) -f -r)8^. c {t)d,(r) = W(f) 

2 JQDi VnK*) Jbd* 

( 10 ) 
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W(r) = -U in (i)-[ -V<m(r)^(r)+ f $,(r-r)0„lMr) ds{r) 

2 JoDi Cfr H r J JdDj 

( 11 ) 

The solution of the above system of integral equations yields values of rp $c 
and d n ip 9C on the dielectric interfaces which can then be used to calculate 
the scattered field at any point r in the exterior domain using Green’s second 
identity: 

Mt) = E j gDi (-G(f - r)d„Mr) + (12) 

The above system of integral equations is of Fredholm first kind and hence not 
well-posed. Nevertheless, the logarithmic kernel allows for accurate numerical 
inversion at the lower incident frequencies and aspect ratios considered in 
the present study. In current work, we are recasting the system into one of 
Fredholm second kind to simulate higher frequencies and aspect ratios. 


3.2 Green’s Functions 


The field at the point r = (x,y) produced by a radiating y-periodic monopole 
of the equation (3) located at the point f = (x,y) is given by the following 
Fourier form: 


C( r-f;*:*) 


_ 1 _ 

47r 


£ 


m=—oo 


%/— 




(13) 


where A m = fc 2 — m 2 , with > 0 = 0, we take 

cj.{yiT^} < o). These inequalities guarantee proper growth and radiation 
conditions in the x-direction. This expression for the Green’s function is 
valid when the parameter values are such that the denominator in (13) is 
nonzero for every m. Since (13) is a slowly converging series when |x - x| < 
1, Ewald’s method 6 is used to re-write the Green’s function in a form that 
is rapidly converging. Via a Taylor expansion in the remaining integral of 
Ewald’s formula, we obtain a representation of the Green’s function entirely 
in terms of special functions that can be easily evaluated numerically: 12 
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+ e-<*-*>^erfc ( e ^ - ^)] (14) 

In (14), E n (z) is the exponential integral function 1 of degree n, R^-(x- 
x) 2 + (y-y + 2n m) 2 and E is an arbitrary real number. The choice E — l/\k\ 
provides an appropriate balance between the singular and regular parts of the 
Ewald representation 12 . In the interior of each rod, we use the free space 
radiating Green’s function of the Helmholtz equation: 

$ j (r-r) = i^ 1) (fc j |r-f|) (15) 

where Hq 1 ^ = Jq + iYo is the zero-order Hankel function of the first kind and 

kj = ky/ij. 

4 Numerical Solution of the Boundary Integral Equations 

For a lattice period consisting of J circular rods, a computational mesh is 
required on the union dD = (Jjfs=i °f dielectric interfaces. We utilize 3- 
node quadratic boundary elements and each rod interface is partitioned into 
M identical boundary elements. We denote the unknowns ift 9C , O n tp ac on the 
mth element of the jth rod by Within each boundary element, a 

local coordinate s (s G [-1,1]) is introduced and the rod surface, the unknown 
fields and the incident wave in (9)-(ll) are all approximated in terms of three 
nodal values via the following quadratic interpolations: 

*4.(«> = Er^m)rK(s) (16) 

tfiW = ELiWJtNAs), d n 4i,(s) = Er.,(S„Vi)r Nr(s) (17) 

The nodal values of the unknowns are given by (V^L)r = V’m( r— 2), (SnV’iJr = 
^nV£»( r ~ 2) and the quadratic basis is: • 

{Ni(s),N 2 (s),N 3 (s)} = ~{ s(l - s),2(l -s)(l + s),s(l + s)} (18) 

To satisfy C 1 inter-element continuity, we move the middle node of each el¬ 
ement to an appropriate point off the circle, the coordinates of which can 
be readily calculated. 12 We substitute (16)-(17) into the integral equations 
(9)-(10), label each mesh node with a global index p, and hence obtain a dis¬ 
cretized system of integral equations at each mesh node r p . Our procedure 
for the numerical solution of the discretized equations consists of evaluating 
all element integrals, assembling a linear algebraic system for the unknowns 
V'p, d n \p p , and solving the resulting linear system. 
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The element integrals involve our two Green’s functions G and $j. In 
evaluating G, it is most efficient to use the Fourier form (13) when |x — £| > 6 
and the modified Ewald form (14) when \x - x\ < <5 where we find that an 
appropriate choice is S = 0.2. The Hankel function in is evaluated us¬ 
ing standard series expansions. 1 Element integrals are computed accurately 
using a standard low order Gaussian quadrature formula or a low order mod¬ 
ified Gaussian quadrature formula with logarithmic weighting function 10 , as 
appropriate. In our computations, all element integrals are evaluated using 
Q =s 8 quadrature points. For more detailed information on our numerical 
integration procedures, the reader is referred to our recent paper. 12 

There are a total of SMJ unknowns which are the real and imaginary 
parts of V'p and d n ipp at each mesh node r p . Once all element integrals have 
been computed, a linear algebraic system can be assembled by passing from 
the local element indices j, m, r to the global mesh index p. Since the matrix 
in the linear algebraic system is dense, we employ a direct solution technique: 
Gaussian elimination with partial pivoting. All computations were performed 
in C ++ on a 233MHz PowerPC G3 processor. 

5 Results 

We simulate experiments that are being performed by H. Everitt and his group 
at Duke University and compute scattering properties for the lattice depicted 
in Figure 1 for the case of normally incident transverse magnetic (TM) waves. 
We compute the transmission amplitude \T(k)\ as defined in Figure 1 for 
reduced frequency k < 1 in which case there is exactly one transmitted and 
one reflected wave mode. All computations assume a rod radius to lattice 
constant ratio of 0.1524 and we consider rods in air with dielectric contrast 
of e = 12.0 for the case of lossless and lossy dielectrics. 

Figure 1 depicts a Fabry-Perot structure that involves two mirrors sep¬ 
arated by a cavity. It is well known that, due to the presence of the cavity, 
transmission through the structure in the bandgap frequency region is possible 
close to resonant bandgap frequencies 8 , as is shown in Figure 4. As the total 
number of rows J increases, the mirrors become more reflecting and the peaks 
in Figure 4 become sharper. At each resonance, a field tends to be maintained 
within the cavity which tends to a localized defect mode as J —> oo. 

We first analyze the lossless case (Im(ej) = 0), starting with the trans¬ 
mission properties of the mirrors themselves, by setting d = 0 and considering 
a square lattice with J rows. The transmission amplitude T is shown versus 
frequency k for a lattice of 4 rows and a lattice of 6 rows (Fig. 2). We have 
found that the number of boundary elements per rod M must be increased 
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0 0.1 OJ OJ 0.4 OJ 0.6 0.7 
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Figure 2. Photonic crystal square lattice with J rows (d = 0). TVansmiusion plots are shown 
for the cases of J = 4 (top) and J - 6 (bottom). Computations are performed with M - 16 
boundary elements per row. 


as k increases and M = 16 was sufficient for the range of frequencies that 
are shown. Similar transmission graphs have been calculated numerically by 
other researchers (e.g. Bell et al. 3 ,Qiu et al. 7 , Sigalas et al. 9 ). 

We see the origin of two bandgaps, roughly one between k = 0.30 and 
k = 0.47 and the other between k = 0.55 and k = 0.68. As we increase J, 
these intervals of k tend to bandgaps of the (doubly periodic) infinite lattice 
in which there is zero transmission. We use the term “bandgap" only in 
reference to transmission of normally incident waves, and not in the sense of 
a complete bandgap which would completely block transmission for all angles 
of incidence. Comparing our results with the atlas of complete bandgaps in 
Joannopoulos et al. 6 , we conclude that our bandgap on the left contains a 
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Figure 3. Fabry-Perot configuration (d > 0). Peak transmission frequencies are shown as 
functions of Fabry-Perot cavity length. The unshaded area corresponds to the first bandgap 
region. 


complete bandgap while the one on the right does not. For our application, 
we locate (approximately) the left bandgap of the infinite lattice by increasing 
J until the frequency interval in which T < 0.01 is uniquely determined to 
two decimal places. Using this criterion, we find that the left bandgap is the 
interval (0.29,0.47). In the context of the finite lattice, we refer to this interval 
as the first bandgap. We also refer to the complement of the bandgaps for 
k > 0 as the set of frequency bands. 

For the photonic crystal with J rows (d = 0), we have observed that the 
transmission amplitude has J — 1 peaks in each band (excluding the peak at 
k — 0). We can trace the location of these peaks through the first band and 
bandgap for increasing values of the cavity length d. An example for the case 
J = 4 is shown in Figure 3. For the range of frequencies shown, accurate nu¬ 
merical results could be obtained with M = 10 elements per row. The results 
clearly indicate that, as d increases, peaks move towards decreasing values of 
k. Each peak that leaves the second band and enters the first bandgap region 
turns into a Fabry-Perot resonance that traverses the bandgap as d is fur¬ 
ther increased, eventually exiting the bandgap and entering the lowest band. 
We have performed similar calculations for J = 8 in the first bandgap (Fig. 
5(top)). The phenomenon of the transmission peaks traversing the bandgap 
persists. 

The horizontal lines in Figure 3 correspond to transmission peaks of each 
of the mirrors (J/2 rows). Based on our observations, there should be J /2 — 1 
horizontal lines in each band. Our explanation of this phenomenon is that 




0.25 0.29 0.33 0.37 0.41 0.46 0.49 

Reduced Frequency (k-fr/c) 


Figure 4. Fabry- Perot configuration (d > 0). Two Fabry-Perot resonances are shown in the 
first bandgap (unshaded area). 


there is total transmission through each mirror at the horizontal line fre¬ 
quencies and, therefore, total transmission occurs through the two mirrors 
successively without any dependence on the cavity length d. 

The sharpness of a transmission peak at bandgap frequencies is of interest 
since, in the case of a laser cavity, it is related to the efficiency of the laser (Fig. 
4). The sharpness or “quality” of the transmission peaks is measured by the 
Q-value which is defined as Q — k/Ak where k is the resonant frequency and 
AA: = |Jk 2 - *i| where = T 2 (k 2 ) = 1/2 (Fig. 5(bottom)). We observe 

a local maximum in each Q-value as the resonance traverses the bandgap. 
New resonances that enter the bandgap peak at successively higher Q-values 
as they traverse the bandgap. When we plot Q-values versus frequency (Fig. 
6), we observe that the maximum Q-values are located at resonant bandgap 
frequencies which are away from the endpoints of the bandgap. 

The curves of Figure 5 (top) are close approximations of the defect modes 
associated with the case J = oo of perfect mirrors. Indeed, the locations of the 
resonant frequencies changes minimally with increasing J . Thus, in Figures 
5-6, we have labeled the resonances DM1, DM2, DM5 where “DM” stands 
for “defect mode” in correspondence to terminology for the case .7 = oo. Our 
numerical results for the locations of the Fabry-Perot resonant frequencies are 
in excellent agreement with the experimental data (Fig. 7). On the other hand, 
Q-values for these defect modes, as predicted by the lossless model, consis¬ 
tently overestimate their experimental counterparts. As such, it is necessary 
to account for dielectric loss in order to accurately quantify laser efficiency of 
the Fabry-Perot cavity under consideration. 
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Cavity Length (d) 
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Cavity Length (d) 

Figure 5. Fabry-Perot Configuration (d > 0). Fabry-Perot resonant frequencies (l.e. 
bandgap peak transmission frequencies) are shown as functions of the cavity length for 
each resonance (top). For the above resonances, Q values are shown as functions of cav¬ 
ity length (bottom). Dots represent computed values and curves shown are cubic spline 
interpolations, while DM stands for defect mode. 


We introduce a small amount of loss into the dielectric rods by adding 
a small imaginary part Im(£j) to the exterior dielectric permittivity in the 
boundary integral model. Our lossy model has been compared to 2-D scatter¬ 
ing experiments performed by H. Everitt and his group at Duke University. 3 
In these experiments, the loss tangent of all the rods is identical with a value 
of approximately tan <5 = Im(ej)/Re(ej) « 0.002. Some representative com¬ 
parisons between the model and experiments are shown for a Fabry-Perot 
cavity with increasing mirror thickness (Fig. 8). Our simulations are in good 
agreement with the experimentally observed Q-values and capture the lim- 
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0J9 0J3 0J7 0.41 045 

Reduced Frequency (k-ft/c) 

Figure 6. Fabry-Perot Configuration (d > 0). For the resonance* of the previous figure, 
Q values are shown as functions of frequency. Dots represent computed values and curves 
shown are cubic spline interpolations. 


7j0| 



Cavity Size (D) 

Figure 7. Measured Fabry-Perot mode frequencies as a function of cavity length for the 
case of square lattice mirrors with J = 2,3,4 and 5 (circles, squares, triangle and diamonds, 
respectively). The solid lines are predictions of the 2-D boundary integral model while the 
dashed lines are predictions of a 1-D layer model. The integer m is used to label the defect 
modes. 


iting behavior with increasing mirror thickness N — J/2 (Fig. 8(left)). It 
PAn be observed that an accurate prediction of Q-values requires both the 
incorporation of loss in the dielectric rods and a genuinely 2-D model. Our 
model dan predicts the drop in peak transmission amplitude with increasing 
mirror thickness that is observed in the experiment (Fig. 8(right)). The de- 
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Figure 8. (left) Measured quality factors (Q-values) for the Fabry-Perot cavity plotted as a 
function of mirror thickness N — J/2 for defect mode m— 4 and cavity length D — d = 3 
(filled circles). Predictions of the boundary integral model with loss (tan S — 0.002,0.004), 
without loss (tan S = 0.000) as well as the 1-D layer model are shown, (right) Peak trans¬ 
mission of this mode as a function of N compared to the lossy (tan S = 0.002,0.004) and 
lossless (tan S — 0.000) boundary integral model. 



Frequency (GHz) Cavity Size ( D) 

Figure 9. Quality factors (Q-values) for the 2-D Fabry-Perot cavity as a function of: (left) 
frequency and (right) cavity length for all modes with J ~ 8(N = 4) from experimental 
measurements (filled symbols) and the 2-D boundary integral model (open symbols joined 
by dashed lines). The grey areas indicate frequencies outside the photonic bandgap. 


feet modes predicted by our model are plotted as functions of frequency and 
cavity length and compared to the corresponding experimental results. (Fig. 
9). It is evident from Figure 9(right) that the model is effective in associating 
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the experimentally observed Q-values with the various defect modes that can 
be sustained in the Fabry-Perot cavity. 

We believe that our boundary integral model serves as a useful tool for the 
experimentalist in that it can quantify the location and quality of bandgaps 
and resonant modes for lossy dielectrics as function of lattice geometry and 
dielectric material parameters. This modeling approach can lead to the for¬ 
mulation of design rules for optimal performance of 2-D photonic crystal de¬ 
vices. In current work, we are extending our boundary integral model of 2-D 
EM scattering to take full advantage of separability in the exterior Green’s 
function. The resulting “fast” boundary integral model will allow for the 
simulation of EM scattering in larger 2-D lattices with localized defects. 
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